1. Field of the Invention
This invention is generally related to methods and systems relating to nuclear magnetic resonance (NMR) measurements and, more particularly, analysis of NMR data using in part a Monte Carlo sampler (or random sampler) that maybe used in oilfield operations.
2. Background of the Invention
Nuclear magnetic resonance (NMR) has been a known laboratory technique and has become an important tool in formation evaluation. NMR well logging background information can be found, for example, in U.S. Pat. No. 5,023,551 to Kleinberg et al., is incorporated herein by reference in its entirety.
In reviewing aspects of NMR, it is known that NMR relies upon the fact that the nuclei of many chemical elements have angular momentum (“spin”) and a magnetic moment. In an externally applied static magnetic field, the spins of nuclei align themselves along the direction of the static field. This equilibrium situation can be disturbed by a pulse of an oscillating magnetic field (e.g., a RF pulse) that tips the spins away from the static field direction. The angle through which the spins are tipped is given by θ=γB1tp/2, where γ is the gyromagnetic ratio, B1 is the linearly polarized oscillating field strength, and tp is the duration of the pulse. Tipping pulses of ninety and one hundred eighty degrees are most common.
It is noted after tipping, two things occur simultaneously. First, the spins process around the direction of the static field at the Larmor frequency, given by ω0=γB0, where B0 is the strength of the static field and γ is the gyromagnetic ratio. For hydrogen nuclei, γ/2π=4258 Hz/Gauss, so, for example, in a static field of 235 Gauss, the hydrogen spins would process at a frequency of 1 MHz. Second, the spins return to the equilibrium direction according to a decay time, T1, which is known as the spin-lattice relaxation time. Because this spin-lattice relaxation occurs along the equilibrium direction, T1 is also referred to as the longitudinal relaxation time constant.
Also associated with the spin of molecular nuclei is a second relaxation time, T2, called the spin-spin relaxation time. At the end of a ninety-degree tipping pulse, all the spins are pointed in a common direction perpendicular, or transverse, to the static field, and they all process at the Larmor frequency. However, because of small fluctuations in the static field induced by other spins or paramagnetic impurities, the spins process at slightly different frequencies, and the transverse magnetization dephases with a time constant T2, which is also referred to as the transverse relaxation time constant.
A standard technique for measuring T2, both in the laboratory and in well logging, uses a RF pulse sequence known as the CPMG (Carr-Purcell-Meiboom-Gill) sequence. As is well known, after a wait time that precedes each pulse sequence, a ninety-degree pulse tips the spins into the transverse plane and causes the spins to start processing. Then, a one hundred eighty-degree pulse is applied that keeps the spins in the measurement plane, but causes the spins, which are dephasing in the transverse plane, to reverse direction and to refocus. By repeatedly reversing the spins using a series of one hundred eighty degree pulses, a series of “spin echoes” appear. The train of echoes is measured and processed to determine the irreversible dephasing time constant, T2. In well logging applications, the detected spin echoes have been used to extract oilfield parameters such as porosity, pore size distribution, and oil viscosity.
In theory, other laboratory NMR measurements may be applied in well-logging to extract additional information about the oilfield, but in practice, the nature of well-logging and the borehole environment make implementing some laboratory NMR measurements difficult. For example, inversion recovery is a common laboratory technique for measuring T1. In an inversion recovery measurement, a one-hundred eighty degree pulse is applied to a system of spins aligned along the static magnetic field in order to reverse the direction of the spins. The system of spins thus perturbed begins to decay toward their equilibrium direction according to T1. To measure the net magnetization, a ninety-degree pulse is applied to rotate the spins into the transverse plane and so induce a measurable signal. The signal will begin to decay as the spins dephase in the transverse plane, but the initial amplitude of the signal depends on the “recovery time” between the one-hundred eighty degree pulse and the ninety-degree pulse. By repeating this experiment for different recovery times and plotting the initial amplitude of the signal against recovery time, T1 may be determined. While this technique has been successfully used in the laboratory for several years, inversion recovery is very time consuming, and those of ordinary skill in the art recognize that inversion recovery may be unsuitable for well logging applications.
Other inversion algorithms are available for analyzing NMR well-logging data. The earliest methods provided one-dimensional T2 (transverse relaxation time) spectra from single measurement data assuming multi-exponential decays. Examples of these methods include the “Windows Processing” scheme disclosed by Freedman (U.S. Pat. No. 5,291,137) and the “Uniform Penalty” method (Borgia, G. C. Brown, R. J. S. and Fantazzini, P., J. Magn Reson. 132, 65-77, 1998) Subsequently, acquisition schemes were devised comprising multiple measurements with different wait-times. Processing techniques were then introduced to analyze these measurements. One such method is disclosed by Freedman (U.S. Pat. No. 5,486,762).
U.S. Pat. No. 6,462,542 issued to Venkataramanan et al. and U.S. Pat. No. 6,570,382 issued to Hurlimann et al. are examples of other NMR methods developed to measure spin relaxation and diffusion in 1D measurements, 2D measurements, and multidimensional measurements. The measurement data is often analyzed by numerical Laplace inversion algorithm to obtain spectra of relaxation parameters, e.g. T1 and T2 and diffusion constant (D). For example, a 1D experiment (such as a CPMG measurement), T2 spectrum is obtained. For a 2D experiment, a joint spectrum of two parameters (e.g. T1-T2, D-T2) is obtained. Several algorithms have been published for 1D experiments, for example by: (1) S. W. Provencher, CONTIN: A General Purpose Constrained Regularization Program for Inverting Noisy Linear Algebraic and Integral Equations, Comput. Phys. Commun. 27, 229 (1982); (2) G. C. Borgia, R. J. S. Brown, and P. Fantazzini, Uniform-penalty Inversion of Multi-exponential Decay Data, J. Magn. Reson. 132, 65 (1998); and (3) E. J. Fordham, A. Sezginer, and L. D. Hall, Imaging ult-exponential Relaxation in the (y; log t1) plane, with application to clay Itration in rock cores, J. Magn. Reson. Ser. A 113, 139 (1995). However, these algorithms are not easily extended to handle 2D data set due to the huge requirement on computer memory.
U.S. Pat. No. 6,462,542 by Venkataramanan et al., discloses new measurement schemes such as “Diffusion Editing”, in which the NMR data is substantially orthogonalized with regard to relaxation and diffusion attenuation, a processing technique based on a separable response kernel has been disclosed (see Venkataramanan, L., Song, Y-Q., and Hurlimann, M., U.S. Pat. No. 6,462,542). This method does not involve any model for the different fluid responses. Instead, it analyses the data in terms of unbiased spectra of relaxation times and diffusion rates. It is attractive in that it requires no a priori knowledge regarding the fluid properties and in favorable cases provides simple graphical results that are easily interpreted even by non-experts. A potential drawback of the inversion is that its accuracy is in part dependent upon the separability of the response kernels. This can limit the range of its applicability to measurements in which the NMR response is substantially orthogonalized in each of the measurement dimensions, for example, application of the method to multiple CPMG sequences with different inter-echo spacings.
Existing processing techniques also impose non-negativity constraints on the individual spectral amplitudes and typically require selection of at least one regularization (smoothing) parameter. The non-negativity condition, based on obvious physical grounds, renders those processing algorithms inherently non-linear. Although not a problem in principle, this places demands on the stability of the selected optimization procedure and caution must be exercised to ensure acceptable repeatability of inversion results for noisy data. The noise issue is addressed by use of the regularization parameter, which ensures that resulting spectra are smooth. However, selecting an appropriate value for the regularization parameter is not trivial. Despite the considerable body of published work addressing the question of regularization from a theoretical point of view (e.g. see references cited in Borgia, G. C. Brown, R. J. S. and Fantazzini, P., J. Magn Reson. 132, 65-77, (1998) and Venkataramanan, L., Song, Y-Q., and Hurlimann, M., U.S. Pat. No. 6,462,542), in practice regularization remains largely subjective, sometimes based only on the aesthetic appearance of the computed spectra. Regularization is of particular importance in multi-dimensional inversions, since the spectra are generally grossly underdetermined by the data and noise artifacts can easily result. In addition, different regions of the spectra display vastly different sensitivity to the input data. Failure to account for these variations in sensitivity can lead to false or unrealistic peaks in the spectra which can easily be misinterpreted.
The inversion of noisy NMR T2 echo data of a T2 spectrum is also widely recognized as an inherently non-unique process (see R. Parker, Y-Q Song, Assigning uncertainties in the Inversion of NMR Relaxation Data, J. Mag. Res. 174 (2005) 314-324.). An approach to quantifying this uncertainty is to use, for example, Monte Carlo sampling. Measurement noise is well described by an uncorrelated normal distribution. When combined with the non-negativity constraint on T2 spectral values, this can lead to spectral values following a non-negative normal distribution. G. Rodriguez-Yam discloses samplers for truncated normal distributions of which non-negative normal samples are a subset, however their algorithms are grossly inefficient for the covariance matrices present in MNR T2 spectral inversion (see G. Rodriguez-Yam et al., Efficient Gibbs' sampling of truncated multivariate normal with application to constrained linear regression, Technical Report, Colorado State University 2004). The reason for this and why it is not a practical method is that they are based on Gibbs' samplers that update the spectral estimate just one T2 component at a time. When all of the spectral elements are fixed but one, that one has little room for change without violating the noise constraints on the data. This means that each spectral sample can only be slightly different from the preceding sample, indicating a high degree of statistical correlation and thus being an inappropriate solution due to the very slow convergence. Thus, what are needed are methods or systems that address all of the above noted problems and among other things, improve convergence as well as open the door for the inversion of 2D NMR spectra.
Accordingly, there continues to be a general need for improved NMR measurements and, in particular for the oil and gas exploration industries, improved NMR methods that can be used to extract information about rock samples and be used in well-logging applications.